Non-adiabacity and large flucutations in a many particle Landau Zener problem 
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We consider the behavior of an interacting many particle system under slow external driving - a 
many body generalization of the Landau-Zener paradigm. We find that a conspiracy of interactions 
and driving leads to physics profoundly different from that of the single particle limit: for practically 
all values of the driving rate the particle distributions in Hilbert space are very broad, a phenomenon 
caused by a strong amplification of quantum fluctuations in the driving process. These fluctuations 
are 'non-adiabatic' in that even at very slow driving it is exceedingly difficult to push the center 
of the distribution towards the limit of full ground state occupancy. We obtain these results by 
a number of complementary theoretical approaches, including diagrammatic perturbation theory, 
semiclassical analysis, and exact diagonalization. 

PACS numbers: 03.65. Sq, 03.65.Yz, 05.45.Mt 



I. INTRODUCTION 

In many physical contexts, one is met with quantum 
systems that are subject to 'slow' time dependent exter- 
nal driving. The most elementary prototype in this cate- 
gory, the Landau-Zener (LZ) system [HH], contains just 
two coupled levels driven linearly in time. In this system, 
the initially occupied instantaneous ground state level of 
its Hamiltonian 

(»-«)■ w 

stays occupied in the infinite future with a probability 

P^l-e-"^. (2) 

where g is the coupling strength, A defines the driving 
rate and exp^irg^/X) is the so-called Landau-Zener pa- 
rameter. The approach P 1 is manifestation of the 
quantum adiabatic theorem, i.e. the statement that suf- 
ficiently slow driving keeps a quantum system in its adi- 
abatic ground state. 

Many quantum single particle systems can be effec- 
tively described in terms of the Landau-Zener setup, or 
one of its multi-dimensional generalizations |31 [H [51 [B] . 
The reason is that the driven approach of pairs of in- 
stantaneous eigenstates will generate 'avoided crossings' 
which can be represented by Hamiltonians such as in 
(jlj. The cumulative statistics of these crossings then 
describes the behavior of the system in the course of the 
driving process. In particular, the system will remain 
in its ground state if only the latter is sufficiently well 
separated from the first excited states. 

But how do many particle systems behave under driv- 
ing? Given the exponential abundance of energy levels 
in interacting systems (or the fact that superimposed on 



the groundstate we often have an continuum of low ly- 
ing, 'soft' excitations) reference to the adiabatic theorem 
will not be sufficient to understand the consequences of 
driving. In some instances, it is possible to reduce the 
problem to one of studying the statistics of linear (oscil- 
lator) excitations superimposed on an invariant ground 
state. A system of this type has been studied in work 
by Yurovsky, Ben-Reuven and Juliene [7j (see also |B.i), 
with the principal observation that the driving process 
generates a number of x excitations, where 

X = exping^X) (3) 

coincides with the Landau-Zener parameter. 

In general, however, many particle systems cannot be 
linearized. The ramifications of non-linearities in a driven 
context have been studied in Ref. [SJITU], within the sim- 
plifying framework of a low-dimensional system (specif- 
ically, a Bose-condensed system which can be described 
in terms of a single complex condensate amplitude.) De- 
scribing its dynamics in terms of a nonlinear Schrodinger 
equation, Ref. |10j observed examples of rather interest- 
ing behavior, including situations where the system does 
not remain in the ground state even in the fully adia- 
batic limit. (For this, however, authors had to consider 
systems of bosons with attractive interactions. In prac- 
tice, such systems tend to collapse.) The problem is also 
particular in that it has a low dimensional Hilbert space. 

Refs. |111 112] considered an interacting driven system 
which is generic in that it shows the two principal char- 
acteristics of interacting quantum systems: a high di- 
mensional Hilbert space and nonlinearity. In one of its 
representations, the system describes a large assembly 
of degenerate fermions which may pair-combine into a 
bosonic level as their energy gets pushed up (a cartoon 
of a fermion/boson conversion process as realized in, say, 
a time-dependent BEC/BCS crossover.) Quite strikingly. 
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it turned out that this system resists getting close to its 
adiabatic ground state, i.e. a state where all particles 
eventually have become bosonic; even at very slow driv- 
ing an 0(1) fraction of particles remains in energetically 
high-lying sectors of Hilbert space. In fact, the meth- 
ods employed in Ref. [TT] did not enable us to get close 
to the 'adiabatic regime' of the model. The approach 
to adiabacity was discussed in Ref. [12^ by mapping the 
evolution of the system to an effective Hamiltonian dy- 
namical system. In this way, as we show in this paper, 
the exponential dependence Q gives way to a power law. 



nb 
N 



A, 



controlling the ground state occupancy for a wide range 
of initial conditions. Here nt, e is the number 

of particles in the bosonic ground state, and N the to- 
tal number of particles. Formally, the reluctance of the 
system to approach its ground state may be understood 
in terms of dynamical instabilities (to be discussed 
in some detail further below.) Physically, it reflects the 
general inertia of interacting systems to adjust to (time 
dependent) environmental changes. 

In this paper, we will discuss an interaction phe- 
nomenon which is no less remarkable: the slow dynami- 
cal evolution of the system goes along with an exception- 
ally strong buildup of quantum fluctuations: although we 
are dealing with a system that should behave 'semiclas- 
sically', on account of the largeness of its number of par- 
ticles, N, the distributions of particles in Hilbert space 
turn out to be very broad. This phenomenon can be 
attributed to an amplification of quantum fluctuations: 
classically, the driving process changes the relative en- 
ergy of two many particle systems. At some point, the 
originally lower (stable) system becomes higher in en- 
ergy (unstable). In a strict classical sense, however, the 
ground state of the elevated system remains stationary. 
It takes the action of weak {0{N^^)) quantum fluctua- 
tions to destabilize this state and initiate the evolution 
towards energetically more favorable sectors of Hilbert 
space. The nonlinearity of that evolution (interactions) 
then leads to an ampliflcation of the initial fluctuations, 
up to a point where they become of 0(1). Specifically, 
we find that 

• in all but the extreme adiabatic limit N — nt, ^ N 
the ground state occupation nf, shows massive fluc- 
tuations, var(nfc)/(nb)^ = 0(1). At intermedi- 
ate mean occupancy (nh) ~ N/2, the distribution 
P(nb) extends over almost the full interval [0, A^]. 

• At fast rates ni, N, the probability distribution 
is exponential: P{ni,) oc exp[— ^^/(nf,)]. At slower 
rates, it becomes even broader and covers Hilbert 
space in a manner for which we have no analytic 
expressions. At yet slower rates, upon approaching 
the adiabatic limit, the distribution gets 'squeezed' 
into the boundary region ni,/N 1. It then as- 



sumes a universal Gumbel form [T3] , the latter fre- 
quently appearing in the context of extreme value 
statistics. 

In this paper, we will set the stage for the discussion of 
fluctuations by establishing contact between our earlier 
'quantum' approach to the problem - the latter adjusted 
to regimes of fast and moderate driving - and a slow driv- 
ing semiclassical formulation. Our semiclassical approach 
is different from that Ref. [T^] in that it microscopically 
connects to the early quantum stages of the evolution, a 
matching procedure necessary to explore the role of fluc- 
tuations. (In fact, we also doubt validity of the results 
for the mean conversion rates derived in that earlier ref- 
erence: for all but zero initial occupancy of the boson 
state, a power law in the driving parameter is predicted 
that we believe is incorrect.) 

The rest of the paper is organized as follows: in section 
[lT| we will introduce the spin variant of the model, and 
discuss the mapping to its other representations. We will 
then analyse the system, in a manner that is structured 
according to the driving rate: in section HI we discuss 



the regime of moderately high driving rates, where the 
system can be effectively linearized and admits a full an- 
alytic solution. In section IV we go beyond the linear 
regime and consider driving rates 1 < A^^ < O(lnA^) in 
terms of effective rate equations. Although these equa- 
tions become uncontrolled for values A"'^ ~ O(lniV), 
they are good enough to signal the system's reluctance 
to approach the adiabatic limit. In section |V] we will 
formulate the semiclassical approach to the model, and 
the so-called truncated Wigner approximation. This ap- 
proximation becomes highly accurate at sufficiently large 
A^ > 10, and for all values of the driving rate. In sec- 
tion IVII we use the method of adiabatic invariants to 
explore the semiclassical theo ry at very slow driving, 
A^^ = 0{lnN). In section [VIl| the quality of the results 
obtained in this way will be checked by comparison to 
direct numerical solutions of the Schrodinger equation of 
our problem (which is feasible for particle numbers up to 
A^ = 0(10^ )) and to the simulations of the semiclassical 
approximation (the latter extensible to much larger A^.) 
We will also compare to previous work in the literature. 
In section IVIIII we discuss a number of ramiflcations of 
our problem relating to previous theoretical and experi- 
mental work. We conclude in section HXl 



II. THE MODEL 

In this section we will introduce the theoretical model 
discussed in much of the rest of the paper - a high- 
dimensional generalization of the spin-boson model. We 
will also introduce a number of less abstract equivalent 
representations, some resonant with concrete experimen- 
tal activity. 
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A. Definition of the model 

Consider an SU(2) spin of value S = N/2 ^ 1, cou- 
pled to a time varying magnetic field of strength (—At) in 
z-direction. This will be the first quantum system par- 
ticipating in the driving process. Its partner system is a 
single bosonic mode at energy {—Xt). We couple these 
two compounds by declaring that the creation of a boson 
goes along with a lowering of the spin by one. The total 
system is then described by the Hamiltonian 

H = -Xtb^b + XtS' + ^{b^S- + bS+), (4) 



N 



where g / \/N defines the coupling strength and — 
Sx ± iSy. The Hamiltonian Q obeys the conservation 
law [i?, 5"= + b^b] = 0, showing that the total Hilbert 
space dimensionality of the problem is 25" -I- 1 = iV -f 1 . 
The linear growth of the dimension of Hilbert space in 
N is of course not representative for 'generic' interacting 
systems (dimensionality exponential in N.) An increase 
in the dimensionality of the problem can be effected by 
symmetry breaking, e.g., by replacement of the large spin 
by an assembly of N non-degenerate spin 1 /2 compounds 
(cf. discussion in section HC below.) While this gener- 
alization will make the problem largely un-tractable, we 
believe it to have little qualitative effect, as long as the 
band-splitting is smaller than the hybridization strength 
with the boson mode. 

Below it will be useful to think of 6 as a transverse 
magnetic field acting on the spin. The above conser- 
vation then implies that the transverse field strength is 
proportional to the deviation of the spin off total polar- 
ization Sz = S. This feedback mechanism of the spin 
precession into the field strength encapsulates the effect 
of interactions in the spin variant of our model. However, 
before proceeding, let us introduce a few alternate rep- 
resentations which make the interpretation of the model 
as one of interacting particles more transparent (cf. Fig. 



1. The variant dominantly discussed in earlier work 
describes the hybridization of energetically degen- 
erate pairs of spinful fermions with a boson mode 
- cf. the Hamiltonian ([t]) below. This may 
be viewed as a dispersionless approximation of 
the fermion-boson conversion processes realized in 
BEC-BCS crossover experiments in fermionic con- 
densates [HITS], cf. Fig. [2]b). (Although the 
initial full occupancy of the flat fermion band as- 
sumed in AG may not be adequate to the descrip- 
tion of the experimental situation [TB], see section 
Vm] below.) 



spin 1/2, the system maps onto a time-dependent 
variant of the Dicke model [17]. In this form it is 
relevant to the description of super radiance phe- 
nomena in molecular magnetism |18j . and to cavity 
QED with many two level systems [T^ . 

3. In a somewhat less obvious representation, the 
model describes the conversion of pairs of bosons 
into dimers. In this incarnation it is of relevance 
to recent experimental work of the JILA group. In 
these experiments, identical [20] or different [2T] 
species of atoms undergo sweep through a Feshbach 
resonance to form di-atomic molecules. (On the 
level of effective classical equations of motions, this 
correspondence was noted in earlier references |1 2] . 
Below, we will establish the connection within a 
fully quantum mechanical setting.) 

For the convenience of interested readers, the equiva- 
lence between these different incarnations of the model 
is established in subsection |II C| below. While the 'spin- 



boson' formulation (4]) does not relate to concrete phys- 
ical systems in an obvious way, we find it ideally suited 
to our theoretical analysis and much of our later discus- 
sion will be formulated in this language. However, it is 
straightforward to transcribe all conclusions to the con- 
text of the other models. 



B. Formulation of the problem 

The problem we will address in much of the paper is 
formulated as follows. Start the dynamics in the dis- 
tant past in the adiabatic ground state of the problem: 
spin fully polarized ^ S ~ N/2 and zero bosons 
nb{t) = {b'^{t)b{t)) = 0. The goal is to find the num- 
ber of produced bosons at large positive times. 



rib = lini nb{t) = lim {b^{t)b{t)). 



(5) 



The adiabatic limit is reached when Ub ^ N or, equiv- 
alently, = —S. Alternatively, we can speak of the 
representation in terms of bosonic atoms and bosonic 
molecules, given by Eqs. (l8| and (10), where this ini- 



2. Identifying the empty (doubly occupied) configura- 
tions of the fermion levels with the two states of a 



tial condition implies the absence of atoms in the distant 
past, and the goal is to find the number of produced 
atoms at large positive times. 

As was anticipated in the introduction, the operator 
b'^b exhibits strong quantum fluctuations. It will, thus, be 
of interest to explore the distribution P{n) defined by the 
moments {{b^b)'^). Also, there may be physical situations 
where it is not appropriate to start the process in the 
spin-up polarized state. We have already mentioned the 
case of dilute fermion gases [T5] where a near south polar 
configuration may be a more appropriate starting point. 
Similarly, a downward polarized state Sz = —S (in an 
initially downward pointing field) will model an atom — > 
molecule conversion process with a purely atomic starting 
configuration. 



4 



In the following, we will describe different stages of the 
conversion process, where the organizing principle will be 
the conversion efficiency, which in turn depends on the 
speed of the sweeping process. For later reference, the 
different parameter regimes, along with the theoretical 
methods we will use to describe them, are summarized 
in Fig.g 



quantum rcgunc 



generic rcgnne 



1. 2. 3. 4. 5. 




adiabatie limit 



FIG. 1: Different theoretical approaches to the problem. The 
Hilbert space of the problem contains a number of different 
sectors: a (quantum fluctuation) dominated region where only 
a few bosons have been created nt = 0(1) and the spin re- 
mains nearly polarized, a 'linear region', where the number 
of bosons may be large but is small enough to justify lin- 
earization of spin operators, 1 ^ ^ A^, a 'generic' region, 
rib = 0{N), and the adiabatie limit nb — > A'". Theoretical 
approaches we will use to approach these regimes include, 
1.) solution of the linearized Schrodinger equation (cf. sec- 
2.) Keldysh perturbation theory based on a self 

3.) semi- 



tion 



nil 



IV I 



consistent RPA approximation (cf. section 
classical approach based on the app roximate conservation of 
adiabatie invariants (cf. section [Vl] ), 4.) an extended semi- 
classical scheme, where semiclassical methods are employed to 
propagate the full quantum distribution beyond the quantum 
regime (cf. sections [v| and VII I, and 5.) numerical solution 



of the Schrodinger equation for moderately large A'^. 



C. Alternate model representations 

The equivalence to a system of boson-coupled degener- 
ate two level systems - the time dependent Dicke model - 



is best seen starting from the latter. Consider the Hamil- 
tonian 



N N 

H = -xt b^b+'^Y.-t + J^i: {b^ + b^t) , (6) 

i—l ^ 2—1 



where a^^y''^ are Pauli matrices acting in pseudospin 
space. The fact that the spin operators appear in to- 
tally symmetric combinations, J^i <t^'^'^ suggests intro- 
ducing the total spin operator S^'^'^ = \ a^'V^^. The 
operators 5'^'*''^ define an SU(2)-algebra. Denoting the 

maximum weight space 'all fermions occupied' by \S), 
we observe that has maximal eigenvalue S = N/2, 
S^\S) = S\S). This shows that the operators S^^-^'^ act 
in a Hilbert space of dimension 25-1-1 = -I- 1. The 
global spin representation of the Hamiltonian ([6| is but 
the starting Hamiltonian The Hamiltonian Eq. ^ 
is represented on Fig. [2] a). As mentioned in the Intro- 
duction, a natural experimental realization of Eq. ^ is 
a cavity QED with N two- level systems. 

The fermi-bose crossover model describes a system of 
2N degenerate spinful fermion states. Two fermions oc- 
cupying such a state may combine to form a boson pop- 
ulating one ('condensate') bosonic state. Starting from 
an infinitely low initial value, the energy of the fermions 
gets pushed up linearly in time (cf. Fig. [2] b) ). This 
setup is described by the Hamiltonian 



H = -Xt 



Xt 



N 

J2 (alt''*! +4 



N 



^ (b^fliia,! + 6a|^aM , (7) 



N 



where A is the driving rate, g/v N the coupling amplitude 
of the conversion process, and a^.n and h are fermion and 
boson annihilation operators. For convenience, we have 
split the time dependent energy symmetrically between 
bosons and fermions. (By a gauge transformation any 
distribution of energies with fixed difference 2Xt between 
a boson and a pair of fermions can be realized.) 



Noting that the interaction couples only two (| f, |) 
and |0)) of four possible (| t,i), I T), I i) and |0)) oc- 
cupation states of a spinful fermion level, we may intro- 
duce a pseudospin state to discriminate between these 
two states, | Tii) (o) and |0) <-> (i). The pseudospin 



representation of the Hamiltonian is but (|6|. 

We may think of the Hamiltonian (|7]) as a cartoon 
version of more realistic models. For example, it can 
be interpreted as a non-dispersive (flat-band) approxi- 
mation to a system of free fermions coupled to the lowest 
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FIG. 2: The different incarnations of the modeh a) time dependent Dicke model: an assembly of degenerate two level systems 
coupled to a boson mode, b) flat band approximation (top) to a system of spinful fermions pair-binding into bosonic 'molecules', 
c) system of atoms pair-converting into molecules, d) spin subject to a magnetic fleld and coupled to a bosonic level. 



mode of an electromagnetic field in a cavity QED setup. 
Alternately, it may be regarded as a cartoon version of 
a BCS/BEC crossover system wherein a band of (lat- 
tice) fermions is initially filled, the band dispersion is 
neglected, and only the lowest condensate state of the 
bosons ('molecules') is kept (Fig. 2]b) inset). While the 
neglect of dispersion effects may be acceptable [T6], the 
assumption of full band occupancy might not be met in 
experiment. (Ref. [T5] argues that the full band assump- 
tion should be modeled by starting the dynamics in in- 
completely 'polarized' (pseudo)spin states, and we refer 
the reader to this work for a discussion of this case.) For 
completeness, let us mention that the problem of the dy- 
namic BCS-BEC crossover has been addressed in many 
publications and using different approaches, beginning 
with t22j and continuing with Refs. [IllISl US [Ml 123 . It 
is probably fair to say that a concensus on whether this 
problem can indeed be reduced to the flat-band model 
with a single Bose mode Eq. ([t]) has not yet been reached 
in the literature. 

Finally, consider a system of atoms sweeping through a 
Feshbach resonance whereafter di-atomic molecules form 
the energetically preferred state (cf. Ref. [20] where the 
conversion processes has been realized in condensates of 
®^Rb or atoms.) This situation is described by the 
Hamiltonian 



H = \td)d - 



\t 



c^c - 



(ic^c^). 



(8) 



This Hamiltonian describes the conversion of atoms to 
molecules, where the former/latter are created by the 



bosonic operators c^/d\ and the symmetry [H,d}d + 
2c^c] = reflects the conservation of the total number 
of particles. The time-independent version of this model 
was extensively analyzed in Ref. [26 . 

The (near) equivalence to the previous model repre- 
sentations is seen by computing matrix elements M^^n = 
{N -n+l,2n-2\H\N -n,2n) , vfhere the state \N-n,2n) 
contains (N — n) d-molecules and 2n c-atoms. Evaluating 
the matrix element, we find 



M 



AT., 



v/(^v^^7rn>(n^^^V2)- 



Now, let \N/2 — n.n)' be the state with n bosons and 
spin polarization 5*2 = N/2 — n. Computing the matrix 
elements M'j^ .^ = {N/2 - n + 1, n - l\H\N/2 - n, n) of 
the Hamiltonian Q, we find 



M'j,^^^^^{N-n + lW 



(9) 



The matrix elements M and M' are identical, up to cor- 
rections of 0{n~^) which vanish in the limit of large par- 
ticle number occupancies. This consideration shows that 
the conversion of atoms to molecules can be discussed 
within the framework of the spin Hamiltonian above, 
where the number of molecules is counted in terms of 
the spin quantum number. 

The correspondence between models can be made per- 
fect by considering heteromolecular conversion processes 
(cf. Ref. 121] for a realization in an ^^Rb — Rb gas.) 
Straightforward generalization of (Is]) obtains the Hamil- 
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model 


system A 


system B 


spin-boson 

time dependent Dicke model 
fermi-bose conversion 
atom-molecule conversion 


spin S = N/2 

N two-level systems 

flat band of 2N fermion states 

di-bosonic molecular state 


boson mode 
boson mode 
boson mode 
boson mode 



TABLE I: survey of the equivalent representations of the model 



Ionian 



\t, 



H = \td}d-'-^{c\ci+clc2) + 



where cj^j creates atoms of species 1/2 and d'^ the (1- 
2) heteromolecule. It is not difficult to verify that the 
matrix elements of this Hamiltonian exactly coincide with 



(d^ciC2 + dcl4), (10) 



those of the spin-boson Hamiltonian ([4]). 

In this incarnation our problem has been analyzed by 
a number of authors, beginning from Ref. [7] (see also 
the more recent Ref. [27^ ) . It is also in this form that the 
problem is the closest to existing experiments [201 HI] • 

Eqs. 0, and ^ define four equivalent defi- 

nitions of the model. For the convenience of the reader, 
the correspondence between the systems involved in these 
definitions is summarized in table HI 



III. LINEARIZED DIABATIC REGIME 

We first consider moderately fast driving rates where 
only a small fraction of the particles undergoes a con- 
version process. In the language of the spin model this 
means that at large times the spin will still be stuck close 
to the north polar regions, lim(^oo('S'— Sz{t))/S <C 1, or 
rib/N ^ 1. In this regime, the curvature of the SU(2) 
spin manifold is not yet felt and the problem can be ef- 
fectively linearized. 

This is best done in the Holstein-PrimakofF represen- 
tation of spin operators 



5- =ct(25-ctc)l/^ 
S+ = {2S-c^cy^^c, 
= S-c^c, 



(11) 



where the algebra of bosonic operators, c and c) de- 
fines the Holstein-PrimakofF bosons. Substitution of 
the linearized approximation S~ = (25*) ^/^c^, 5*+ — 
(25) ^/^c, S"" = 5 - c^c into ^ generates the quadratic 
Hamiltonian 



iJnp = -\t {h^h+ c^c) + g{b^c^ + cb). 



(12) 



To compute the number of bosons created in the lin- 
earized evolution, we switch to the Heisenberg represen- 



tation h{t) 
value 



be 



iHt 



and consider the expectation 



nb{t) = {b\t)h{t)), 



taken with respect to the vacuum state, (...) = (0| . . . |0) , 
with 6|0) = c|0) = 0. (Alternatively, one can consider 



a model of bosonic atoms and molecules, in the regime 
where most of the p artic les are molecules. In the Hamil- 
tonians ([s]) and Eq. (10) one may then replace d \/~N , 
to arrive precisely at Eq. (12 1, an appropriate relabeling 
of operators understood. In Ref. ^ a model equivalent 
to the linearization above was introduced, and the mean 
conversion rates of the linearized dynamics were com- 
puted (see also Ref. J8j). 

The computation of n?, = limt^oo nb{t) is reviewed in 
Appendix with the principal result 

Tlh — X — \. (13) 



Eq. (13) states that the number of converted particles 
grows linearly in the LZ parameter. Unlike with the stan- 
dard (two-level) LZ problem, the many particle system 
remains parametrically detached from adiabacity, even 
for large values of the LZ parameter. However, the lin- 
earization underlying the derivation limits the validity of 
the result to n\j <^ N . 

The techniques used in Appendix |A] actually yield the 
entire (quantum) distribution of the number of produced 
particles, 



X 



-n/x 



(14) 



This distribution is very broad. Its width ((n^) — 
(n)^)^/^ ~ X is of the same order as its mean. This is 
remarkable inasmuch as our system is not subject to any 
sources of fluctuations, other than its intrinsic quantum 
fluctuations. We are met with the unusual situation that 
quantum fluctuations are of the same order as mean val- 
ues, in spite of the 'semiclassical' largeness of the latter, 
(rt) ~ a; 3> 1. Similar phenomena have been discussed in 
the context of dynamical sweeps through quantum phase 
transition points, cf. e.g. Ref. [28] . 
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where 



IV. THE KINETIC EQUATION 



D" oF -FoD' 



A. Keldysh formalism 

An alternative way of analyzing the problem starts 
from the boson- fermion representation (TJ) and analyzes 
the corresponding kinetic equation - theory strand b) 
in Fig. [2] To formulate this approach, we employ 
the Keldysh formalism and introduce the appropriate 
Keldysh Greens functions (cf. Ref. [53] for a review of 
the formalism.) Readers not interested in the formalism 
are invited to proceed directly to the discussion of the 
main result of this section, the rate equation (27 1. 



Within the Keldysh formalism, the fermion subsystem 
is described by the matrix Greens function 



act. a' a' 



where a, a' — 1,2 index the two-component Keldysh 
space (to which the matrix structure in the equation 
refers), and a = is a container index comprising 

the fermion site and spin index. G^, G^, and G^ stand, 
respectively, for the advanced, retarded and Keldysh 
Greens function. The latter is defined by 

where the symbol 'o' defines a temporal convolution, 
iAoB){t,t') = f dt" A{t,t")B{t'\t') 



and f = I —2nf defines the distribution of the fermions. 
hi what follows, it will be convenient to represent the 
Greens functions of the theory in the Wigner representa- 
tion. 



G^(t,c.) 



dt G^ 



t 



t 



(15) 



and similarly for all other Green functions. In the ab- 
sence of interactions, g = 0, a straightforward calcula- 
tion obtains the Wigner representation of the free Greens 
functions as 

GLAr,^) - J^,. (16) 

GL.'(t,co) = -2Ti5,„.fo{r,\t/2)8{Lo-Xt/2), 

where /o(t, Ai/2) = /o(t) describes the distribution of 
the full fermion band, /o(t) = 1 — 2no(T) — —1. In a 
similar manner the bosons are described by 

D{t,t') = ~iK{t)bl{t')) {t,t'), (17) 



and F = 1 + 2rn,, rib describing the boson distribution. 
The Wigner representation of the non-interacting boson 
Greens functions is given by 



(18) 



1 



LO + \t 

D^{t,Lj) = ~2TlFo{T,-\t)5{LJ + \t), 



where Fo(t, —At) = F[){t) relates to the vanishing boson 
number of the non-interacting problem as Fq{t) = 1 + 
2no(r) = 1. 

In order to calculate the number of produced bosons at 
the end of the Landau-Zener process, we need to calculate 
the exact bosonic Keldysh Greens function, in terms of 
which 



Ub = lim nb{t) = lim (iG^{t,t) - l) /2 



(19) 



where Ai, = —2lmD^ is the spectral function of the 
bosons. For large positive times, the Greens functions 
become effectively uncoupled, D^{uj) asymptotes to the 
free form ( 18 1 and Ai,{uj) ~ 2TrS{uj + Xt). This means that 



the number of bosons at the end of the process is given 

by 



rib 



lim 



F(r, -At) - 1 



(20) 



B. Perturbation theory 



First, we can try to calculate the bosonic Greens func- 
tion perturbatively, in powers of g (which means, in prac- 
tice, in powers of the dimensionless parameter g^/X). 
This expansion is valid at fast rate A, g^/X ^ 1. The 
lowest order correction to the bosonic propagator is de- 
picted in Fig. [3] diagram a). The calculation of this di- 
agram, although cumbersome, proceeds in a straightfor- 
ward fashion, giving the result 



Hb 



A 



(21) 



Each vertex of the diagram carries a factor of g/vN. 
However, there are N fermions propagating around the 
loop, canceling the factor 1/iV coming from the verices. 
This results in the cancellation of factors of TV in the 
answer. 

In the next order of perturbation theory two more di- 
agrams contribute to the bosonic propagator. These are 
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b) 



including g and A kept fixed. Our starting point is the 
system of kinetic equations (the Wigner transforms of the 
Dyson equations, cf. Ref.f29') for the boson and fermion 
distribution function. 




d) 



e) 



{dr - Xd^) F{t,u;) = I (S^ - (E« o F - i^^ o E-^)) (t, w), 

(24) 

(d, + /(r, e) = * (fi^ - (a^ o / - / o a^)) (r, e), 

(25) 



FIG. 3: The lowest order corrections to the bosonic propa- 
gator. 



shown on Fig. |3| dia grams c) and d). This leads to the 
second order answer 



rib 



A 



1 

N 



(22) 



Note that at iV = 1, this indeed coincides with the Taylor 
expansion of the exact Landau-Zener formula 



rib 



1 



(23) 



Notice also that at > 1 , the number of produced par- 
ticles is bigger than at = 1. This is natural, as with 
larger N one can produce more bosons. 



ate 



Finally, we note that the calculation leading to Eq. ( 22 ) 
closely parallels that employed in Ref. 24J to calcu 
the molecule production in a Feshbach resonance exper- 
iment perturbatively. 



The first of these equations is obtained by evaluat- 
ing the Wigner transform of the general kinetic equa- 
tion of a boson system [29^ {-idr + [H, ])F = - 



(E^ oF-FoT."^) for the Hamiltonian H = -Xt. 
second equation is its fcrmionic analog. 



The 



In the 'coUision integrals' of ((24|, E^^-^^^ (cr^-^"^) are 



the components of the bosonic/fcrmionic self energies. It 
is not difficult to see that the dominant contributions to 
these self energies are given by the diagrams Fig. [s] a) 
and e), external legs truncated. Unlike with the lowest 
order perturbation theory discussed above, the fermionic 
and bosonic propagators appearing in these diagrams are 
to be understood as dressed propagators, containing self 
energy insertions by themselves. (Technically speaking, 
this means that we will evaluate the self energies in a self 
consistent RPA approximation.) Self energy diagrams 
with crossing interaction lines, such as Fig. |3]b), carry 
factors of relative to the self consistent RPA dia- 
grams. (See, however, the discussion in the end of this 
section.) 



C. Self consistent RPA and kinetic equation 

Our aim now is to calculate the number of produced 
bosons in the large N limit, with all other parameters 



The derivation of the self energies for the boson- 
fermion interaction follows standard procedures |29j and 
we will not repeat it here. As a result we obtain the 
equations 



where we suppressed the explicit r-dependencies for no- tational clarity and A — 



-2ImG^ and Ab = 



-2Imi:>^ 
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are the spectral functions of the fermions and the bo- 
son, respectively. Notice the factor of N^^ multiplying 
the collision integral of the fermions. The absence of 
this factor in the boson collision integral signals that the 
boson interacts with all N fermions simultaneously. In 
contrast, each fermion interacts only with a single boson, 
which means that the fermionic self energies carry the 
uncompensated squared coupling constant g'^/N. 

Eq. (26 1 defines a particle number conserving sys- 
tem. Consider the quantity N{t) = ni,{T) + Nnf{T). 
Using Eq. (19 1, and its fermionic analog nf{T) — 
— \ (/ ||/(t, e)A(r, e) — l), we obtain that the number 
of particles varies in self consistent RPA as 



dN 1 
d7=2^^' 



de 



Integrating the first (second) equation (26 1 over 
/ '^Ab{T,uj) (/ ^A{T,e)) and using that Xd^A^ = 
drAf,, we find that the two integrals cancel out, and 
drN = 0. 

This is about as far as we will get without further ap- 
proximations. Below we will argue that the validity of the 
RPA approximation is limited by the condition Ub <^ N. 
In the limit nt/N 0, the self energies cr^/^ broad- 
ening the energy dependence of the spectral functions 
become vanishingly small. (The detailed dependence of 
the broadening depends on the value of time, r.) Relying 



on this limiting behaviour, we will approximate the spec- 
tral functions by 74(e) ~ 27r(5(e) by (5-functions. This is a 
bold approximation inasmuch as the pairs of (5-functions 
in (26 1 enforce the resonance condition t — 0. However, 
it is this 'resonant' time window, r ~ where the distri- 
bution functions entering the integral kernel are expected 
to vary strongest. This means that the value of the col- 
lision integral may well be sensitive to the detailed tem- 
poral profile of the distribution functions around t — 0. 
Ultimately, we will need to compare to the result of other 
methods to verify the validity of the ^-function approxi- 
mation. 

We substitute the approximations Ai,{t,uj) ~ 2TTd{uj + 
At) and A{T,e) ~ 27r(5(e - At/2) into ([26]), use ^{.(t) = 



l {J^FiT,u;)Ab{T,uj)-l) ~ (F(t,-At) - l)/2, and 
differentiate w.r.t. time to obtain 



drUi, = —drUf 



A 



5{T)[n} + nb{2nf-l)]. (27) 



The meaning of this equation becomes transparent upon 
rewriting the combination of distribution functions on its 
r.h.s as (1 -|- ni,)n'j — ni,{l — nj)^. This factor weighs the 
probability that two fermions convert into bosons and 
back. To solve the equation, we take advantage of the 
conservation law 

nt + Nnf=N. (28) 
It is then straightforward to obtain the solution 



nbir) 



-e(r) ^N{N + 4) +2^ + N+ .yN{N + 4) 



(29) 



The singular behaviour of n;, at t = is a remnant of the 
above (5-function approximation. Incidentally, we note 
that the Taylor series expansion of with respect to 
gives 



Hb = lim nb{T) 



A 
6 





4) 











8iV + 
67V2 



A 



.(30) 



The first two terms indeed match the direct perturba- 
tive expansion given by Eq. ( 22 1 , confirming that this 



technique does take into account all the diagrams up to 
the second order. (In the context of perturbation theory, 
the (5-function approximation of the spectral functions 
amounts to the lowest order in g Born approximation to 
the self energy operator.) 



On the other hand, taking the limit iV — s- oo gives 
ivfe^-1^ 

rib = 



2e^ +N 



(31) 



This is the main result of the kinetic equation approach. 
For e^a^/^ < iV we find 



nfc = e 



1. 



(32) 



This matches the large N limit given by Eq. ( [T3| . 

The kinetic equation produces the boson number (31 1 
which is monotonously increasing with decreasing A, un- 
til it reaches its limiting small A value of N/2. The rea- 
son for this behaviour is that at large N , the r.h.s. of 
the kinetic equation approaches zero at ny = 1/2. In 
other words, the rate of fermions converting into bosons 
is matched by the rate of bosons converting into fermions. 
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However, we need to remember that at ni, = N/2 we 
are well outside the regime of validity of the approxi- 
mated kinetic equation. In fact, we are outside the limit 
of the RPA as such. This is because at rib ~ dia- 
grams that are nominally small in become sizeable, 
on account of the large value of the distribution func- 
tion njj N . For example, it is straightforward to verify 
that its crossing interaction lines make diagram Fig. [3] 
b) small in iV~^, but this smallness gets counteracted 
by n^-dependent propagators in the center regions of the 
diagram. The bottom line is that the RPA is controlled 
only in the limit ni,/N — > 0. Our comparison with numer- 
ics below will show that it works remarkably well beyond 
its nominal limits, i.e. the profile (31 1 represents a good 
approximation even at ni,/N = 0(1). However, in order 
to advance into the adiabatic regime, rif, — s- N, we need 
to employ different methods. 



V. SEMICLASSICAL ANALYSIS. 

Our goal now is to understand the behavior of n?, for 
arbitrary driving rates A. We will use that for large N ^ 
1 our problem approaches a well defined classical limit 
where both spin and the bosonic fields can be treated 
as classical objects. The classical equations of motion, 
obtained by replacing operators in the Hamiltonian Q 
by c-numbers and commutators by Poisson brackets, read 



ib = 
S^, — 



Sy — 



-Xtb 



N 



s- 



-i^ibS^ 



b*s-), 



-i-^{b* ^b)S^- XtSy 
{b + b*)S, + \tSa,. 



(33) 



Within the strictly classical setting, these equations have 
to be solved for the initial conditions: b{to) = 0, Sz{to) = 
N/2, Sx{to) = Sy{to) = 0, where to — oo is the 
initial time. It is straightforward to verify that the 
solution corresponding to these initial conditions reads 
{b{t),S,{t)) = {b{to),S,ito)) = (0,A^o), i.e. the initial 
configuration of zero bosons represents a classically sta- 
tionary (if unstable) solution. 

The prediction that no bosons will be generated, is 
clearly wrong. However, the purely classical treatment 
has yet another, and related, drawback: for given initial 
conditions the uniqueness of the classical solution does 
not permit the buildup of fluctuations. This is in contra- 
diction with our earlier findings that the distribution of 
rib is wide at least at fast and intermediate rates. 

It is thus obvious that an meaningful description of 
the classical limit must account for the presence of quan- 
tum fluctuations. These fluctuations will lead to an ini- 
tial destabiliziation of the configuration {b{t), Sz{t)) — 
{0,No). The nonlinearity of the classical Hamiltonian, 



then amplifies the fluctuations and drives the system 
away from its initial state. The expansion of dynamics in 
quantum fluctuations for bosons was analyzed in detail 
in Ref . [30' . Using a path integral construction similar to 
the Keldysh technique it was shown that to leading order 
in quantum fluctuations the classical (Gross-Pitaevskii) 
equations of motion do not change. However, the initial 
conditions become randomly distributed with the weight 
given by the Wigner transform of the density matrix. 
This so called Truncated Wigner Approximation (TWA) 
originated in quantum optics |31l 132] and has later been 
adopted to cold atom systems [33] • By by now the TWA 
has become a standard tool in describing the dynamics 
of cold atoms in regimes where classical Gross-Pitaevskii 
equations do not suffice (see Ref. [3?| for review). 

The TWA is most conveniently adopted to our model 
by using the Schwinger boson representation of spins. 
The small parameter of this expansion is 1/5 where S 
oo describes the classical spin. In Appendix |B] we give 
some details of this derivation. The net result is that the 



classical equations of motion (33 1 need to be supplied 



by stochastic initial conditions distributed according to 
the Wigner function. For large spin (pointing along z- 
direction) and vacuum state of bosons, representing the 
initial conditions in our problem, this Wigner function is 
approximately given by 



ttS 



exp[-2rto] cxp --^ S{Sz-S), 

(34) 

where Sj_ = ^ S'^ + Sy and uq is the initial boson num- 
ber. The fact that S± and uq are not exactly equal to 
zero is the result of vacuum quantum fluctuations. In this 
work we are primarily interested in the number of bosons 
created during the dynamical process, and its fluctua- 
tions. The connection between moments rip, 1 = 1,2,... 
and the variables b and b^ is established by the Wigner 
transform identities summarized in Appendix [B] For in- 
stance. 



no 



b*b- 1/2, 



{b*bf-b*b. 



(35) 



(The correction term 6*6 appearing in the expression for 
the second moment ensures that zero point fluctuations 
do not affect the expectation value of no and Ug, which 
should be zero in the vacuum state. In the classical limit 
rio ^ 1, these terms become inessential.) 

By way of an illustration, let us outline how the TWA 
reproduces the results of the quantum calculation in 
the linearizable regime. For quadratic Hamiltonians the 
TWA is exact [33 [31] and we must be able to repro- 
duce all the results of Sec. jlllj Upon leaving the linear 
regime, the number of bosons n{t) = (6^^(i)6(i)) has be- 
come large and semiclassical methods are expected to 
work with high accuracy. In other words, the exactness 
of the TWA to the linear regime entails its applicability 
to the full parameter space of our problem. 

In the linear regime, we may assume approximate con- 
stancy of 5"^ : 5"^ « 5' = N/2. Eq. ^ then reduces 
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to 



ib = —Xtb + gs 
is^ = —gb + Xts^ , 



(36) 



where we introduced s~ = S~ /yN. These equations are 
identical to an exact reformulation of the Schrodinger 
equation (cf. Eq. (Al I in Appendix [A| , a consequence of 
the linearity of the problem. In Appendix |B 3| we show 
that the solution of these equations, with account for 
quantum fluctuations in the initial data, indeed repro- 
duces the results of the full quantum calculation. In 
particular we show there that rih = x — \ and (n^) = 
— 3a; + 1 (here as before, we reserve the notation rib 
for the number of bosons at the end of the process, as 
in Eq. Both results perfectly agree with the exact 

distribution ([l4|. 

At slow values of driving, x ^ N the nonlinear equa- 
tions (33) no longer afford an exact solution. In the next 



section we will discuss an approximate solution scheme, 
based on the method of adiabatic invariants. In Sec. IVIII 
we will analyze in detail nh{X) and its distribution in 
various regimes using both exact and TWA simulations. 



VI. DEEP ADIABATIC LIMIT 

A. Adiabatic Hamiltonian 

In this section, we will apply the concept of classical 
adiabatic invariance to the adiabatic limit of the driving 
process. Considering a fixed initial value no ~ 1 (con- 
sistent with Eq. (34)), we will be sloppy about initial 
conditions. This is good enough to obtain reliable re- 
sults for the mean value of produced bosons, (nf,), but 



won't suffice to obtain the statistics of the driving pro- 
cess. The latter will be explored in section |VII| below 
within the more sophisticated TWA scheme. 

The concept of the adiabatic evolution of classical in- 
variants [35J was previously applied to the semiclassical 
limit of a non-linear Schrodinger equation |10j and later 
[T^ to the classical equations of our model. However, 
(for all but very small, compared with TV, initial oc- 
cupancy TiQ ~ 1) the latter reference predicts a power 
law A^/^ for the conversion rate which we can exclude 
on the basis of direct numerics, TWA, and our analyti- 
cal calculation below. We therefore believe that the ac- 
tual technical calculation of Ref. [I^ is flawed. However, 
the conceptual basis for the emergence of power laws (as 
opposed to exponentials such as in the original Landau 
Zener problem) remains valid: in driven classical dynam- 
ical systems, certain invariants of the autonomous limit 
of the evolution (viz. the action picked up upon traversal 
of periodic trajectories) become weakly time dependent. 
While in most regions of phase space this time depen- 
dence is very (exponentially in the driving rate) weak, it 
can become strong (algebraic) in the vicinity of singular 



points. Such points reflect the presence of nonlinearities 
in the Hamiltonian, which in our problem are due to in- 
teractions. When the driving parameter sweeps across 
a singular point, the topology of trajectories in phase 
space (and thence the action picked up along these tra- 
jectories) may change profoundly. It is the dynamics of 
these processes which we will investigate in the following. 

Specifically, we will consider the system in the follow- 
ing limit. 



Af > 1, 



A 



9^ 



«1, 



1 



<4logA^<l, (37) 



where a > is an arbitrary positive number. 

We consider the Hamiltonian (|4| in the limit of c- 
number valued operators. Using a number-phase de- 



composition of the boson field, b 



and a po- 



lar representation of spin variables, Sz = Scos9, Sx = 
5 sin 61 cos Sy — 5 sin sin ^, it assumes the form 



N 



-Xtn+Xt— cos 9+gy Nnsin (6) cos(^— (/j). 



The two angles ^ and (p can be combined into one angle 
(p = ^ — if + TT (where tt is added for later convenience) . 
The angle 9 and the boson number n are not independent, 
but are related via the conservation law 



n = — (1 — cos 



(38) 



This allows to trade 9 for n. It is then convenient to 
rescale 



Nn, 



(39) 



so that the new n varies from to 1, with the last 
value taken when all particles are bosons. Similarly, the 
rescaled initial condition reads 



no 



1 

N' 



It is also advantageous to replace 



H NH. 



(40) 



Throughout this section, we will work with the rescaled 
variables as described here, restoring the original vari- 
ables at the end. 

Finally, without any loss of generality we can set 5=1 
since it can be scaled out by the appropriate rescaling of 
t. We arrive at the Hamiltonian 



H 



-jn — 2n\/l — ncos(0), 7 = 2Xt. 



(41) 



Here (f> and n play the role of the coordinate and momen- 
tum, canonically conjugate to each other. The equations 
of motion of these variables read as 



dtn 



dnH. 



(42) 
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These equations of motion need to be supplied with ini- 
tial conditions. In this section, where the emphasis is 
on the calculation of the 'typcial' number of produced 
bosons, we will supply the dynamical system (42) with a 
fixed value of the initial n. We take this to be the typical 
value no ~ 1 of the (rescaled) distribution ( 34 ) . (As 
we will see below in the slow limit n(t — s- oo) has only 
logarithmic sensitivity to ng and precise form of the ini- 
tial value is not important, as long as we are interested in 
estimating mean conversion rates. For the same reason 
the initial fluctuations in the angle 9, which we ignored, 
are also unimportant.) 

We begin our analysis of the initial value problem by 
plotting lines of constant energy H = const, at differ- 
ent fixed values of 7. Examples are shown in Figs. [4] [s] 
|6] and [Tj The lines shown in these figures represent 




FIG. 4: Trajectories of If at 7 = —3. Here the horizontal 
axis is n, and the vertical axis is (p- 





FIG. 6: Trajectories of iif at 7 = 1 




FIG. 7: Trajectories of Ji" at 7 = 4. 



trajectories of the Hamiltonian for 7 fixed. In our prob- 
lem, 7 changes in time. However, at small A, 7 changes 
slowly. Thus the system will follow one of the trajecto- 
ries for some time before it will slowly drift to a different 
trajectory due to the slowly changing 7. 



B. Fixed Points 

An important role in our analysis is played by the fixed 
points of the Hamiltonian system [101 HZ] • This is where 
the trajectory consists of a single point (7 still kept fixed.) 
They are found by solving 



dH 



0, 



dH 

dn 



0. 



(43) 



There are two families of solutions (n, (f)) of these equa- 
tions. The first exists for 7 > —2 and has (p — 0, 
{n,(p) = (711(7), 0), where 



FIG. 5: Trajectories of J? at 7 = —1.3 



'^1(7) 



12-72 + 7^12 + 72 



18 



(44) 
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This solution appears at ni(2) = at 7 = —2. As 7 
is increased above —2, ^1(7) starts increasing until it 
asymptotically reaches lim^_>oo ?t-i(7) = 1- 

The other solution has (f> — t:, and requires 7 < 2. It 
is given by 



"-2(7) 



12-7^-7^12 + 72 
18 



(45) 



This solution appears at 712(2) = 0. As 7 is de- 
creased below 2, 7x2(7) starts increasing and asymptotes 
as lim^^_oo "-2(7) = 1- 

We thus arrive at the following picture. Initially the 
number of bosons is given by = no ~ 1 /N, at undeter- 
mined 4>. The coordinates n and (j) evolve according to 
the approximate Hamiltonian H ft! —7/1, where we used 
that at large negative times 7 is very large to keep only 
the leading order contribution in this scale. This means 
that the motion is initially given by 



1 



= -lt. 



(46) 



However, as 7 is increasing, the fixed point ^1(7), given 
by (44 1, appears at the critical value 7 — —2. As 7 keeps 
increasing the point moves towards larger values of n. At 
a certain moment of time, the trajectory (n,(j))(t), will 
drastically change: it will no longer perform extended 
propagatation in i/i-directions but start winding around 
the fixed point Eq. ( 45 ) . Eventually, at some large posi- 



tive 7, the trajectory will stop winding around the fixed 
point and merge into a new extended trajectory, which 
at large positive 7 will take the form 



rib, = -jt. 



(47) 



71b is the final number of bosons we are trying to calculate. 

The fixed point 711(7) ^i^s an intimate relationship with 
the instantaneous ground state of the quantum Hamilto- 
nian Eq. Q. This is discussed in more details in Ap- 
pendix ic) 



C. Adiabatic invariants 

In the classical mechanics of slowly parameter depen- 
dent Hamiltonians there exist approximately conserved 
'adiabatic invariants' |35j . Generally, adiabatic invari- 
ants assume the form of action integrals ~ J pdq, where 
the integral is over one full revolution of the system's 
motion. In the present context, the invariant is given by 



1 

2n 



d(j) n. 



(48) 



To further elucidate the meaning of the adiabatic invari- 
ant in our problem, let us go back to the angular variable 
e according to Eqs. ([38]), ((39|. We find 



/ = 



1 

47r 



'(1 



COS! 



(49) 



This is nothing but the area of a part of the surface of a 
unit radius sphere bounded by the trajectory 0{(p), where 
9 and cj) are thought of as spherical angles, divided by 47r. 
This geometric definition of the adiabatic invariant im- 
plies an ambiguity which needs to be resolved. A closed 
trajectory on the surface of a unit sphere separates the 
surface into two domains, and one needs to pick the area 
of one of them. 

We define the adiabatic invariant as the area (divided 
by 47r) of the domain on the surface of the unit sphere 
bounded by the trajectory such that the trajectory, when 
the domain is viewed from above, encircles the domain in 
the counter-clockwise direction. We note that this defini- 
tion does not always conform to the algebraic expression 
Eq. (48 1. However, it results in the adiabatic invariant 



continuous under smooth deformations of the trajectory, 
while Eq. ( 48 1 can be discontinuous under those trans- 
formations. 



For trajectories such as the initial trajectory Eq. (46 1 



or the final trajectory Eq. (47 1 the adiabatic invariant is 
very easy to calculate. Initially the trajectory encloses 
the north pole of the sphere in the counter-clockwise di- 
rection, to result in 



initial 



7(7 ~f —00) 



1 

2^ 



d(j)n{-j — > —00) 



1 

N' 



At the end of the process Eq. (47) implies that the trajec 



tory now encloses the south pole in the counter-clockwise 
direction. The adiabatic invariant is now given by 



i final 



= 7(7 00) = 1 — 



1 

2^ 



d(j)n{"/ 



In the limit A ^ 0, / is conserved. 



'final 



1 - "-6 ^ 



1 

N' 



00) = 1 — rift. 

(50) 

(51) 



This implies that at the end oft the process all particles 
become bosons, Tit, « 1 (with 1/N factor coming from a 
quantum uncertainty irrelevant in the large A'^ limit). 

As the rate A is increased, the adiabatic invariant starts 
changing with time. In the next section, we will review 
the general theory of adiabatic action changes, and then 
apply it to our particular problem at hand. Readers not 
interested in details of the adiabatic dynamical evolution 
of the system may just note the final result of our anal- 
ysis, Eq. (l63l. 



D. 



The theory of approximate conservation of the 
adiabatic invariants 



The theory of approximate conservation of adiabatic 
invariants is a well developed subject, described in detail 
in Ref. f35^. One of its main results is that the change 
of an adiabatic invariant during some time dependent 
process is usually exponentially suppressed at slow rate. 



A/. 



exp 



const \ 



(52) 
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where A is the rate of the process. In fact, one may inter- 
pret the standard Landau-Zener transition probabihty in 
terms of this behavior. 

However, the theory also states that if a singularity 
develops at a certain time during the evolution, Eq. (52) 



breaks down and gets replaced, typically, by a power law. 
We are going to see that this is indeed the case in our 
problem. 

For the convenience of the reader, we briefly summa- 
rize the essentials of the theory of adiabatic invariants, as 
discussed in chapters 49 and 50 of Ref. [35 . Consider a 
Hamiltonian H = H{<f>, n, 7) depending on a pair of con- 
jugate variables (</>, n), and a slowly varying parameter 7. 
For starters, assume 7 to be constant. Since H describes 
a system with only one degree of freedom, we have inte- 
grability, and it is convenient to transform to canonical, 
or 'action-angle' variables. The action variable is defined 
as 

where the integral is over a closed trajectory of the sys- 
tem's evolution. Such curves are specified by a value of 
the conserved energy, E, and, in our case, by the value of 
7. We thus have a functional relation / = I{E,'j). Now, 
consider the so-called abbreviated action 



S{^,E,j) 



d^'n{E,(^\^), 



where the integral extends only over a certain segment of 
the trajectory. The relation / = I{E, 7) may be inverted 
to express 5* as 

5(0,i?(/,7),7) = 5(0,/,7) 

as a function of the coordinate and the action variable. 
The abbreviated action is the generator of a canonical 
transformation (0, n) —f (/, w) from the old coordinates 
to the (conserved) action, and a conjugate angular vari- 
able. Specifically, we have the relation n = d^S{(t>, 
and define the angular variable as w = djS{(j), I 
The equations of motion in the new variables read 



dtl = - 



dH' 



dtw = 



dw 
dH' 



where H' = H'(I,w,"f) is the Hamiltonian expressed in 
terms of new variables. For an autonomous system (7 = 
const.), H'{I,w,-f) = H{I,j) = E{I,-f) is but the old 
Hamiltonian expressed in terms of the action variable. 
In this case, the action is conserved, and the angle varies 
as u; = tdjE. During each period of the motion, the 
action changes by an amount AS = 27r/. The relation 
w — diS then implies a change of w by 2tt. Thus, 



w{t) = t 



27r 

Y' 



For a non-autonomous system, classical mechanics states 
that 

H'{I,w,^) - F(/,7) + ^ = F(/,7) + ^da, 

where S on the r.h.s. must be expressed as a function of 
{I,w) after the parameter differentiation. In other words. 



H'{I, w, 7) = H{1, 7) + A(/, w, 7)dt7, 



where 



A(/,u;,7) = 



55(7, 0,7) 



97 



(53) 



0— 0(/,'u;,7) 



This then means, that the action changes, and it does so 
according to 



dtl 



dA{I,w,j) 
dw 



da- 



rn 



The change of the action over a finite interval of time 
then follows as 



A/ = 
where 



-2A 



dw 9A(/, w, 7) 



dw 



(55) 



2tt 

Y 



is the characteristic frequency of the trajectory. For time 
independent 7, the integral is over one full interval of a w- 
periodic integrand and vanishes; it is temporal variations 
in 7 that give it its finite value. The relations above 
express the evolution of the action entirely in terms of 
quantities that do not contain explicit time dependence. 
One may argue that AI evaluate to Eq. (52), unless 



the time dependence of and R contains singularities. 
This will turn out to be the case in our problem. 

Figs. |4] [5] [6] and [7] show examples of trajectories that 
enter the computation of action integrals for our dynam- 
ical system. (Each of these trajectories corresponds to 
a particular value of / while w parameterizes the trajec- 
tory.) In computing the above integrals, we will then be 
met with the following scenario: consider a trajectory, 
close to the left boundary of phase space and at an ini- 
tial value 7 <C 0. Initially, that trajectory will be 'open', 
i.e. its angular variable will periodically run from to 
2tt, at moderately changing n. To the trajectory we may 
assign an instantaneous value of energy H(n, 0, 7), which 
changes on account of increasing 7. As 7 approaches the 
value —2 from below, we run into a singularity: at a cer- 
tain value 7 = — 2-|-eo;0<eo<Cl, the energy vanishes, 
H{n, (p, —2 -I- eo) = 0. That point marks a drastic change 
in the topology of the trajectory, it changes from open to 
closed. See Fig. [s] for an example of a closed trajectory, 
and a 'critical' trajectory with H = 0. The latter begins 
and ends at n = 0, and has the form of a horizontally 



15 



aligned horseshoe. The passage time through the critical 
trajectory is infinite, and it is logarithmically divergent 
on the closed trajectories for e \ eo (in contrast, it can be 
checked that another horseshoe-shaped trajectory, which 
starts and ends at n = 1 and can be seen on Fig. [7] is not 
critical and its passage takes finite amount of time). The 
action integral receives its dominant contribution from a 
range of e- values above eg. 

To describe this situation algebraically, we expand the 
Hamiltonian Eq. (41 1 for small e, n and (f), where 



7 = -2- 



as 



H 



The singular trajectory is given by _ff = 0, or n = e — 
0^. Its frequency is zero (the period is infinity) and its 
adiabatic invariant is given by 



1 

2^ 



2 3 
62 . 



(56) 



The trajectories inside the singular trajectory have H < 
(one of these is shown in Fig. |5]) . On these closed tra- 
jectories, n is a two valued function of 4>, where the two 

branches n±(0) = ^ (^e ^ (fi'^ ± (^{(j)^ - e)^ + AH^ ^'"^ 

represent the right and left bending arc as cf) varies in 
-\l < </) < \/e-2y/^. The action of these 

curves equals the area enclosed by the trajectories, di- 
vided by 27r, or 



1 

2tj: 

£3/2 



'4 
3 



H 

' 72 



e)^-H4iJ ; 



log 



H 



16e2 



- 1 



(57) 



This expression is approximate, valid for small | <C e^. 
The corresponding frequencies are given by 



OJ ■ 



dH 



log(^ 



H 
16e2 



(58) 



These frequencies indeed vanish logarithmically as H 
approaches zero, confirming that the H = trajectory 
is critical. 

The mechanical system we need to study starts its evo- 
lution ai t —^ —oo by moving according to Eq. (46 1. Sub- 
sequently, at a time where 7 is slightly larger than —2, the 
system crosses the singular trajectory and starts moving 
around the fixed point Eq. (45 1. It is at this time that 



the adiabatic invariant receives most of its increase. 

We estimate the increase by first calculating S for the 
part of the trajectory represented by n^^cj)) where < 



4'2 - ef + 4H 



In this formula, H = H{I) must be understood as a 
function of /, which can be found by inversion of (57 1. 



We then calculate A according to ( 53 1 



A 



1 f"^ 

2 70 ' 



where (j) has to be substituted hy (p — 4>{I, w, A) after the 
differentiation. We now use this result to compute the 
time variation of the action according to ( [54| ) . Using that 
the w-dependence of A is in (j){w, /, A), and that d^j — 2A, 
we have 



dj ^ ~2\--^ = -\ 
0(p ow 



1 - 



' de 



dw 



Eventually, this quantity has to be integrated over time, 
and to do so, we need information on the time depen- 
dence of 4){I,w,^) = (/)(/, wt, 7). The equations of mo- 
tion tell us 



dt4> = d.nH=-^{e-cp^) 



AH. 



(59) 



We next assume that during the time when the invariant 
/ receives most of its increase, takes values such that 
the r.h.s. of this equation vanishes. The validity of this 
presumption will be checked in the end of the section. 
This implies 



(60) 



Indeed, everything done so far implies \H\ <^ and 
so the right hand of Eq. ( 59 1 approximately vanishes if 



Eq. (60) holds 



Using this stationarity condition, and dw(j) = dt(i 
we obtain 



In this equation H — H{I, e) ^ is implicitly defined 
by (57). This gives 



dH 
Ik 



ar / di 

Ih \dH 



2e 



In I 



H 



a term that is larger than the \J —H contribution to dtl. 
Using (58), we arrive at the estimate 



dtl- 



2Aei/2 



Using that e = 2 -I- 2 At, we integrate this expression over 
e from e = eqi remembering that / — 2tJ^ jiZ-n) when 
e = eo, to find 



^3/2 



37r 



(61) 
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This expression looks formally identical to Eq. ( 56 ) , but 
has a different meaning. While Eq. (56) is the adiabatic 



invariant of the critical trajectory, Eq. (611 signifies that 



this expression remains approximately true even at later 
times when the system no longer follows the critical tra- 
jectory. 

We finally need to find the maximum value e = e* for 
which Eq. (61 ) still holds (at even higher values of e, the 



adiabatic invariant will stop growing and just oscillate 
about its average value). The criterion we will use to 
determine that value reads w = tt, i.e. we demand that 
the system proceeded along half a period of the critical 
trajectory. At larger instances of time (e) the transit 
into the domain of oscillatory motion and no further sys- 
tematic action increase has taken place. In practice, the 
fixation of e* turns out to be somewhat tedious, and we 
have relegated it into appendix [D] As a result we obtain 



-A log TV 



(62) 



Finally, with the help of Eq. (511 this leads to our final 
answer 



nb~iV( 1-— logA^ 



A 



(63) 



where we have re-introduced the coupling strength g and 
returned to the original (unrescaled) rib. Note that in the 
derivation above we relied only on the fact that the un- 
rescaled no ~ 1 due to quantum fluctuations. If we start 
from the state where no is non-zero (but smaller than N) 
due to initial occupancy of the molecular state the result 
( 63 1 remains valid with the only difference that the argu- 
ment of the logarithm should be changed to N / uq. This is 
qualitatively different from the results of Ref. [12 where 
an instant crossover to a ~ A^/"^ power law the moment 
'^o ^ 1 /N was predicted. 



VII. TWA AND COMPARISON WITH 
NUMERICAL SIMULATIONS 

In this section, we will apply a combination of numer- 
ical diagonalization and the TWA to obtain accurate re- 
sults for both, the mean value of n^ and its distribution. 
For moderate values of N, the Schrodinger equations 
for Eq. ([4| can be solved directly. Indeed, the Hilbert 
space corresponding to that Hamiltonian, in the sector 
S^ + b^= N/2, has only iV + 1 states. Thus the Hamil- 
tonian reduces to a {N + l)-dimensional matrix. This 
matrix can be diagonalized at reasonable numerical cost 
up to values N < 10"^ (see Appendix [c] for an example of 
how this procedure can be set up). 

At larger values of A^, we simulate the classical equa- 
tions of motion, with initial drawn from the quantum 



may obtain results for significantly larger values of N. 
In the following we will apply both methods interchange- 
ably. However, before doing so, let us first test the accu- 
racy of TWA by comparison to numerical diagonalization 
at values of N where both methods are applicable. 



A. TWA vs direct diagonalization 

In section [V| we have argued that the exactness of the 
TWA in the linear regime entails its applicability to the 
whole range of driving parameters. To back up this claim, 
let us compare TWA results to those obtained by numer- 
ical diagonalization. In Fig. [sjwe show (n;,)(A) for two 
values of A^ = 25": N = 64 and A^ = 128. The sohd 
and dashed lines represent the exact and TWA solutions, 
respectively. There is no visible difference between them. 
To demonstrate that TWA is not actually exact we show 
in the inset the difference between TWA and numerical 
results multiplied by a factor of 100. Clearly for N = 128 
the accuracy of TWA is better which signals that as N 
increases TWA gives results of increasing precision. 



-Exact- — - TWA 



N=128 

Exact TWA 




Wigner distribution Eq. (34 1 - the TWA. In this way, we 



FIG. 8: Dependence of the scaled number of bosons rib/N as 
a function of the rate A. The two pairs of curves correspond 
to iV = 2S = 64 and iV = 128. The solid lines represent 
the exact solution while the dashed lines are the semiclassical 
TWA. The inset shows the difference between exact and semi- 
classical results (A^ = 64 red solid line and A'^ — 128 orange 
dashed line) multiplied by a factor 100. 

In Fig. [9] we plot th e relative number fluctuations 
Srib/N = VK) - {uby/N as a function of A for A^ = 64. 
As with the dependence nfc(A) the difference between 
TWA and exact results is minuscule and it vanishes fast 
as A^ gets larger. These results make us confident that 
TWA represents is high precision method which becomes 
asymptotically exact in the limit N 00. 
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FIG. 9: Dependence of the scaled variance of the number of 
bosons 5nt/N after the process as a function of the rate A for 
= 64. The two curves represent the exact and the TWA 
results. 



FIG. 10: Number of remaining bosons Ub vs the rate A for 
different values of A^. It is clear that as A'^ increases the depen- 
dence goes from exponential form characteristic for Landau- 
Zene r p rocess to the linear one predicted analytically (see 
Eq. j63|) 



B. Combined application of TWA and direct 
diagonalization 

In the following we will analyse our system by com- 
bined application of the two numerical methods. We will 
not indicate which method was used for which particular 
curve, unless necessary. In all numerical simulations we 
set g = 1. (This can be always achieved by the rescaling 
A A(7^.) Unless stated differently, the notation nb(A) 
refers to the mean value of the number of bosons. 

In Fig. 10 we plot nt,{X) for different values of N. For 
small N the function nh(X) is qualitatively similar to an 
exponential form: rtf,(A)/iV ocl — l/a;=l — exp(— vr/A), 
as with the conventional two-level Landau-Zener prob- 
lem. However, as N increases this exponential behavior 
disappears and the dependence becomes much closer to 
the linear one predicted in Eq. ( |63[ ) . To demonstrate the 
approach to linear scaling, we plot ni,{X) for large values 
of TV up to iV 



10 in Fig. 11 The dashed line is the fit 
to the linear dependence: rib/N « 1 — 0.28Aln A^. Note 
that the prefactor 0.28 is slightly different from the I/tt 
in Eq. ( 63 1 . But given that the analytic result was ob- 



tained using a series of approximations we find this level 
of agreement to be satisfying. 

Fig. [12] shows the results of numerical simulations at 
N = 10^ N = 10"*, N = 10*^ and N = 10^ plotted as 
(1 — nh/N)/logN vs \/{TTg'^). The straight line in the 
plot represents the deep adiabatic limit Eq. ( 63 1 , while 



N=10' 
N^IO" 
N=10' 
N=10' 
1-0.28 /I ln(N) 




1.5 2.0 2.5 3.0 



?tln(N) 



FIG. 11: Same as in Fig. ( 10 1 but for large values of A^. Note 



that the horizontal axis is AlnA'^. 



the second continuous line is the solution of the kinetic 
equation Eq. ([Sl]) at = 10^ only. 

We see that the deep adiabatic relation Eq. ( 63 1 



IS 

confirmed in the limit of slow rates. As to the solution of 
the kinetic equation, while it describes well the relatively 
fast rates, it deviates from the numerics at slower rates 
and eventually breaks down when it saturates at nonzero 
value when A — > 0. 

Next, let us discuss the distribution function P(ni,). In 
the fast limit we have a well defined analytic prediction 
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l-Hb/N 




2 



P{no) = 2exp[— 2no] we can invert Eq. (64) and derive 
the distribution function of ni,: 

P{nb) « 2^C2e-P exp [-2C72A^e-''] , (65) 



Ci / rib 



The function P{nb) above describes the Gumbel distri- 
bution, which often appears in the context of the extreme 
value statistics. 



FIG. 12: Number of remaining bosons rih vs the rate A for 
different values of A''. Here circles represent A'^ = 10'', squares 
iV = 10'', rhombi = 10*^, and triangles A'' = 10*. The 
continuous lines represent our analytic solutions, as discussed 
in the text. 




FIG. 13: Distribution function of nj, for A'^ = 128 and two 
values of the rate: A = 2 and A = 4. The distributions are 
discrete and the lines are the guide to the eye. The dashed 
lines, which are barely distinguishable from the numerically 
computed solid lines, are the analytic prediction ( 14 ) 



given by Eq. p4| While for the slow limit we do not have 
an analytic prediction for the distribution function, one 
can make a good ansatz based on the result Eq. ( 63 1 . 



Note that this result was obtained for a particular choice 
of no = 1. More generally for each fixed initial value of 
no it should read: 



nb{no)/N » 1 



A 



ln(C2iV/no), 



(64) 



where Ci is a constant close to tt and C2 is another 
constant, which in general can depend on A. Find- 
ing this constant is beyond accuracy of our analytic 
derivation. Because we know the distribution of no: 
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FIG. 14: Distribution function of for A' = 128 at 
A = 1, 0.4, 0,2. The solid lines are the result of the semi- 
classical TWA approximation and the dashed lines are the 
distributions obtained by exact diagonalization. The dotted 
grey line (which almost coincides with the solid blue line) is 
the fit to the Gumbel distribution ( |65[ ) with Ci ~ 3.05 and 
C2 « 0.03.) 



We plot the distribution function P{nb) in Fig. 14 
The solid lines are the result of the TWA approxima- 
tion and the dashed lines are the exact result. At slow 
rates the distribution function becomes highly asymmet- 
ric. At A = 0.2 it is very well fitted by the Gumbel 
distribution (65) (grey dotted line) with Ci ~ 3.05 (in- 
stead of tt) and C2 ~ 0.03. The fact that C2 is so small 
probably indicates that this constant depends on A. We 
checked that the fit also works well for N = 256 with 
the same constants Ci and C2 . From Fig [14] it is appar- 
ent that the semiclassical approximation correctly pre- 
dicts the Gumbel-like shape of the distribution function 
at small A but misses the oscillations on top of this shape. 
These oscillations have a minor effect on the expectation 
value of Hb and on its fluctuations. However, they are 
not negligible and they do not vanish if we increase N 
(at least they persist up to A^ = 512). We are not sure 
what the origin of these oscillations is. It is interesting to 
point out that a similar interpolation of the distribution 
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function from the exponential to the Gunibel form was 
found in the completely different context of interference 
between two independent Luttinger liquids [3^ . 



Discussion 



As a result of a straightforward adaption of the discus- 
sion of section |IV[ we find that the number of bosons nf, 
for Eq. Q in the infinite future behaves as 



N 



rib 



(66) 



The result ( 63 1 implies that the thermodynamic limit 
of our system never behaves adiabatically: no matter 
how small the driving, in the limit N —^ oo, the adia- 
batic limit becomes elusive. Such type of behavior has 
been found earlier |37j for more general low-dimensional 
bosonic systems near instabilities such as second order 
phase transitions. The giant quantum fluctuations ac- 
companying this non-adiabatic behavior result from the 
inversion of the classical ground states of the participat- 
ing systems. Zero point fluctuations are required to trig- 
ger a macroscopic inversion of state occupancies and the 
stochastic nature of this initiation process generates vast 
fluctuations at later stages. 

Finally, it is worthwhile putting the results above in a 
more general context: a number of previous works [38l 
IMJ HOI SH SI SI mi have addressed the slow dynamics 
of low-dimensional critical systems. In all these cases it 
has been argued that the number of excitations produced 
upon slow driving scales as a power law of the rate (as 
opposed to the the exponential dependence Eq. ^ for 
the Landau-Zener problem.) The exponent of this power 
law is related to the critical exponents characterizing the 
phase transition [38j |39] under consideration. In low di- 
mensions with their relatively higher density of low en- 
ergy states there is a stronger tendency for smaller pow- 
ers, i.e. less adiabatic regimes. In a sense, our 'infinite- 
range interaction' or 'mean field' system exemplifies these 
structures in the extreme setting of a zero-dimensional 
system. 



VIII. APPLICATION TO OTHER DRIVING 
PROCESSES 

We spent most of this paper discussing driving pro- 
cesses where the bosonic sector of the system was initially 
empty (or there were no a atoms in the language of the 
atom-molecular conversion systems Eqs. ([sf and p^.) 

For these initial conditions, quantum fiuctuations are 
needed to jump start boson (or molecule) production, 
and this is what ultimately leads to the power law ( 63 1 . 



We here briefiy discuss the 'inverse' problem, where ini- 
tially all N particles were bosons (or atoms). This vari- 
ant was previously discussed in Ref. [27]. What makes 
it different, is that no quantum fluctuation are required 
to start the conversion process. Solution of the classi- 
cal equations of motion Eq. (41 1 at inverted driving rate, 
7 — — 2Ai, and with initial condition n{t — > — oo) ~ N 
captures the dynamics of the process and quantum fluc- 
tuation introduce only minor corrections. 



At fast rates, this is equivalent to rib — Ne ^ , A ^ 
TTg^, which is but the standard Landau-Zener predic- 
tion. As the rate slows down, below x = e^^ Z-*' ~ 1 
the Landau-Zener prediction breaks down. (Notice that 
the initial [linear] regime in the direct problem was much 
more robust.) 



Taken at a face value, Eq. (66) predicts that at slow 



driving rates ni, saturates at 1 /2. This cannot be correct 
since Eq. (66 1 must break down at slower rates. A de- 



tailed analysis involving the quasiclassical Hamiltonian 
Eq. (41 ) with the initial condition uq = I, gives 



7if,/A^~ 0.057A, 



(67) 



where the constant of proportionality 0.057 was deter- 
mined numerically from the linear flt (note that unlike 
that in Eq. (63 1, it is A^-independent) . In Fig. 15 we 




FIG. 15: The number of remaining molecules vs. A for the 
inverse problem, where we start from the full molecular con- 
densate and no atoms. Note that the results for A*' = 32 
and A*' = 64 almost identically coincide. The thin red line 
represents the linear fit at slow rates: nt/N ~ 0.057A.) 

show numerical data for the dependence nf,(A) for the 
inverse problem for two values of iV = 32 and N — 64. 
The thin red line represents the linear fit at slow rates: 
rib/N « 0.057A, and the numerical data approaches this 
asymptotic in an TV-independent manner. 

Thus the number of remaining bosons, a measure of de- 
viation of adiabaticity, remains proportional to the rate 
at slow rates, as opposed to the exponential suppression 
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of the standard Landau-Zener process. Even though one 
approaches the adiabatic regime when the Landau-Zener 
parameter is close to 1, the deviation from adiabaticity 
is still significant and much larger than in the standard 
problem. This may have significant implications for the 
experiments such as those reported in Refs. [201111], pro- 
vided the initial state is chosen to be atoms and one 
is interested in the molecule production. We also note 
that, as far as we can see, this point was not discussed in 
Ref. [37]. 

In some applications it may be important [16 to ex- 
plore the sector of the theory where the total number of 
particles is much smaller than N (i.e. much smaller than 
the number of two-level systems in Eq. ([6|, or smaller 
than the number of fermionic pair levels in Eq. (|7|.) 
Once more, it is easiest to analyze the system in its 
spin-iV/2 incarnation Eq. Q. Translated to the lan- 
guage of this Hamiltonian a sparse initial occupation of 
fermionic Hilbert space means that we start the process 
at zero bosons and a spin that is pointing nearly down- 
wards (while the total spin continues to define the di- 
mension of state space, S = N/2.) The particle conser- 
vation law then reads Sz + n — —N/2 + e, where e <^ N 
and < n < e is the number of bosons. This regime 
can be modeled in terms of Holstein-Primakoff bosons. 
However, as the spin continues to move downwards, the 
regime of Holstein-Primakoff linearizability will not be 
left. Introducing 



For small values of the LZ-parameter, 



5+ = ct(25-ctc)i/2 
S- = (25-ctc)i/2c, 
= -S + c^c, 



(68) 



where S = N/2 (the version of Eq. (Ill for the spin 
pointing close to downwards), we find 

H=Xt - c^c) + g {b^c + ch) . 

Here h^h -\- c = e. This problem is equivalent to the 
standard Landau-Zener problem. This proves that in this 
sector the model reduces to the standard problem and 
does not have any novel features. 



IX. CONCLUSIONS 

Let us summarize the findings of our paper. The model 
described by Eq. (as well as its other incarnations, 
given by Eqs. (6]), iQ, and Q together with (lOl), with 
the initial conditions 



lim (b^t)b{t)) = 0, 

results in a total final number of bosons n?,, defined by 
nb= lim (bUtWt)) 

given by the following expressions: 



< N 



1. 



This result can be obtained within the Holstein- 
Primakoff method, or by the large-iV approxima- 
tion within the Keldysh formalism. 



• At intermediate values. 



N 



rib 



iv(e^-l) 



2e- 



N 



This can be found using the kinetic equation ap- 
proach within the Keldysh formalism. 

• And at large values, 



> N 



Ub^Nll- 



-\og{N) 



This result was found using the quasiclassical ap- 
proximation together with the adiabatic invariants 
formalism. 

We observe main distinctions between our results and 
what would have happened had we extrapolated the stan- 
dard Landau-Zener to our problem (which would have 

given Ub — N(l — e ~)-) First, the transition between 
the fast rate to slow rate regime happens at 



A 



logiV' 



the rate which becomes progressively smaller as TV is in- 
creased (unlike the Landau-Zener result A ~ tti?^). Sec- 
ond, if A <C TTg^/logiV, then Ub approaches A'' as a lin- 
ear function of A, much more slowly than the Landau- 
Zener result which would predict an exponentially fast 
approach. Finally, particle occupancies in our system 
show massive quantum fluctuations, comparable in value 
to the mean particle numbers. 
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APPENDIX A: DERIVATION OF EQS. ([I3j) AND 

<P|) 



Consider the Heisenberg equations of motion i dtx - 
\H,x\, X — c,b corresponding to the Hamiltonian (|l2|. 



idtb+Xtb + gc' = 0, 
-idtc^ + Xtc^ + gb=0. 



(Al) 
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The structure of these equations suggests the linear com- 
bination 



b = pbo- 



(A2) 



where cq = c(to), io ^ is some time in the distant past 
when the dynamical evolution is started, and /i, . . . , cr 
are complex valued functions depending on time. Dif- 
ferentiating the first of the two equations in (All once 



more w.r.t. time and substituting the second equation, 
we obtain 



[d^ -iX+ {\tf ~g^]x^Q, a; = c, h. 



(A3) 



These two equations are coupled by the initial conditions 
which follow from the first order ancestor equations ( Al I 



To make these dependencies explicit, we substitute the 
expansion (|A2[) and obtain 



[3j — i\ + (Xt)^ — 5^] K — 0, K — iijV, p, cr, 

K = p,p: K{to) = l, dtK{to) = iXto, 

K = v,a : K(to) — 0, dtnito) — ig. (^^4) 



Notice that 



fj. = p and i> = a, 



(A5) 



where the first line is a consequence of (/^, p) and (v, a) 
obeying differential equations with identical initial condi- 
tions, and the second is enforced by the fact that the time 
dependent operators {c,c\b,b^) obey canonical commu- 
tation relations. 

These equations afford a solution in terms of parabolic 
cylinder functions [8 . All we need to know about these 
functions is that 



2 t^oo TrgVA 



Using rib — (c'c ) 



i^bo){p,co + ybl)) 



(A6) 
^ and 



Eqs. (|A5]), (|A6| we then obtain Eq. p^ . 

The actual distribution of the particle number can be 
obtained by a bit of linear algebra. We begin by con- 
densing Eqs. (A2) into the matrix equation 



V p 



Co 



Co 



p ~v 

— V p 



where we used Eqs. (A5). Denoting the occupation num- 



ber eigenstates of the time dependent operators by |fc, 
i.e. 

c\kj) = Vk\k-l,l), c^\kj) = VkTl\k + l,l), 
b\k,l) = Vl\kJ-l), b^\k,l) = VTTT\k,l + l), 



We now expand the vacuum state (of the operators 
co,6o), 1 0), in the basis {|fc, Z)}: 

|0> -^Cfc,z|fc,/). 

k,l 

To determine the coefficients Ck^i, we use the defining 

relation = co|0) ^ {pc ~ vb^)J2k,iCk,i\k,l)- This 
equation, and its partner 6o|0) = generate the set of 
recursion relations 



co|0) = : pVk+ Ick+i^i+i = Wl + lck,h k,l>0, 

pck,o = 0, fc > 1, 



6o|0) : i/Vfc + IckA = pVr+lck+i,i+i, fc, Z > 0, 

pcoj = 0, 1>1. 

From these equations we deduce c^.q — for fc > and 
Co,; = for I > 0. The first and third equation then 
imply Ckj — for k ^ I. The diagonal terms successively 
descend from co,o according to 



Cfe ; — Ski co.o 



The normalization condition 1 = (0|0) = J2k i gen- 
erates the additional condition 



1 



|co,o| 



1 - 



|co,olVr 



We thus arrive at the expansion 



1^1 ^To 



From this, the moments of the particle number operator 
are straightforwardly obtained as 



(0|(ctc)'"lO) 



1 °° I 

„|2 |„ 



\2k 



1^1 ^0 



2k 



fc" 



Comparison with the formulation of moments in terms 
of the discrete probability distribution P{n) 



Cm = ^ P{n) 

n 

leads to the identification P{n) 



(A7) Eq. ( |A5| , 



P(n) 



or, usmg 



(A8) 



Substitution of Eq. (A6) then leads to Eq. (14). From 



this distribvition, the mean value is readily obtained as 



(131 
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APPENDIX B: DERIVATION OF THE 
TRUNCATED WIGNER APPROXIMATION FOR 
THE SPIN-BOSON MODEL. 

1. Truncated Wigner approximation 

In this appendix, we first review the TWA for a gen- 
eral system of bosons, and then adapt to our specific 
spin-boson problem. Consider a (generally time depen- 
dent) Hamiltonian formulated in terms of normal ordered 
bosonic creation and annihilation operators Ti{'4'\ip,t). 
The hat-notation ijj is used to distinguish operators from 
c-numbers. In the classical limit the operators ^p are sub- 
stituted by complex numbers ip satisfying the Hamilto- 
nian equation of motion (Gross-Pitaevskii equation): 

* dt S^p* 

where the right hand side denotes a variational deriva- 
tive. Classically, these equations have to be supplied by 
definite initial conditions ('0Oj V'o)- Within the truncated 
Wigner approximation - the first order quantum correc- 
tion to the classical picture - the initial data becomes a 
distribution W{ipo, iPq) dipodipQ. The kernel of this distri- 
bution is defined by 



X e 



(Bl) 



where p is the initial density matrix of the system. Ex- 
pectation values of operators f2('0^, i/') at time t are then 
to be calculated as [30] 



m)) 



dijjQdroWo{ibo,ro)^ci{m,r{t)), 



(B2) 

where tp{t) is the solution of the classical Gross-Pitaevskii 
equations of motion with initial condition ipo f^nd 
Uci{ip,ip*) is the Weyl symbol of the operator fi. For 
a normal ordered the latter is defined by 

n,i{ij,r)^ J J dT]di^*nii;-T]/2,i^* + ri*/2) e-'^l'/^, 

(B3) 

where il{ip,^*) is obtained by substitution of operators 
in fl{'4', ip^) as '0 ^ "0 ip''^ ip* . For example for the 
number operator — we get 



0*1/' 



1 

2' 

(B4) 

The same result can be obtained by first symmetrizing 
number operator with respect to "0^ and 



4. - '0''^'0 -|- 4>1p'^ 1 



and then substituting the operators 0^ and ip with the 
numbers ip* and ip (see Ref. [34]) for more details. 



2. Adaption to the spin-boson problem 

Our next task is to generalize the TWA to the Hamil- 
tonian To this end, we employ the Schwinger repre- 
sentation of spins through bosons a and f3: 



■/3t/3 



S- = j3h 



(B6) 



Here, the boson operators a and (3 satisfy the additional 
constraint n = a^a + f3^f3 = 2S. The operator n com- 
mutes with all spin operators which means that the ful- 
filment of the constraint at an arbitrary time will be con- 
served for all later times. The option of a purely bosonic 
representation means that the TWA can readily be gen- 
eralized to the Hamiltonian ^ . Once the TWA has been 
formulated for the Schwinger bosons, we are at liberty to 
switch back to spin variables Eqs. (B6 1. However, at this 



point, the spins have become classical numbers, defined 
in terms of the c-number valued boson variables a and 
p. 



The classical equations of motion are given by Eq. ( 33 1 



The initial density matrix describes a pure spin state, po- 
larized along the z-axis or, in bosonic language, a state 
with 25' a bosons and no (3 bosons. (For the spin pointing 
along — z a and (3 should be interchanged). The corre- 
sponding Wigner function then reads [30) [3l] [45] : 

W{a,a\P,n - 2e-2|"l'-2|/3pL2s(4|a|2), (B7) 

where Lm{x) is the Laguerre's polynomial of order N . At 
large S, the latter is strongly oscillatory and the Wigner 
transform is localized near |ap = 25-1-1/2 (see Ref. [50]!. 
So in this case to a very good accuracy (up to 1/5^) one 
can use 

W{a,a\P,P*) « y2e-2|/3|'j(|Q,|2 _25- 1/2). (B8) 

If we reexpress this distribution function through spins 
then to the same accuracy we will find 



■0' l/) 



(B5) 



where Sj_ = 3^ + Sy. This Wigner function has a trans- 
parent interpretation. If the quantum spin points along 
the z direction, because of the uncertainty principle, 
the transverse spin components still experience quantum 
fluctuations so that 

which is indeed the correct quantum-mechanical result. 
Finally, the distribution of the 6-boson represents the 
Gaussian Wigner function of the vacuum state, 

W{b,b*) = 2exp[-2fe*6] ee W{n) = 2exp[-2n], 

where in the second representation we switched to a rep- 
resentation in terms of the boson number n = h*b. 



23 



3. TWA vs quantum solution of linear regime 

We here discuss how the solution of the equations of 
motion ( 36 1 obtains information equivalent to that of the 



full quantum solution of the linear regime. To this end, 
we notic e that the equations afford a solution as (see 
Eq. (|A2l) 



oo we have = 



bit) = iib{to) 



li*s-{tfi) + y*h{ta), 

X and \v'^\ = a; — 1 
(see Eqs. (A5l and (A6l. If we are interested only in the 



where at t 



Dcr of bosons we do not need to know 



statistics of the num 
the phases of /i and v. Using the Wigner function (Bl) 



to compute the average and the Weyl symbol ( 35 ) for the 



operators rifc, it is then straightforward to obtain 



X — 1, 



(nfc) ^ x{nb{h)) + {x- l)(|s-(io)P) - ^ 
and the second moment 

+ Ax{x - l){nb{to)){\s-{tQ)\^) = 2x^ - 3a; + 1. 
It is easy to check that this result conforms with the 



distribution ( 14) 



APPENDIX C: THE MEANING OF THE FIXED 
POINTS 



Now consider a semiclassical version of this Hamilto- 



nian, Eq. (41 1, which for completeness we reproduce here 



H = —jn + 2ny/l — ncos(0). 

Here both the Hamiltonian and the variable n were 
rescaled according to Eqs. (39 1 an d ([40| to ease compar- 
ison with the language of section |VI[ although we need 
to keep that in mind when comparing it with Eq. ( C2 1 



which was not rescaled in any way. Let us minimize this 
Hamiltonian with respect to n and (j). Its minimum is 
given by the substitution of the fixed point = 0, and 
ni('-f) given by ^1(7) — for 7 < —2 and by Eq. (44), or 

12-72+7^12 + 72 



'^1(7) = 



18 



if 7 > — 2. To find the energy minimum we substitute 
71,1(7) i^^to H to find 



-N 



Emin ^ NminH ^ 0, if 7 < -2 

367-73 + (12 + 72)? 
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if 7> -2.(C3) 



The minimum of the classical Hamiltonian Eq 
together with the eigenvalues of Eq, 
as a function of 7 on Fig 



C3f, 

(C2|, are plotted 
One can see 



[16] for iV = 10 
that the classical minimum closely follows the quantum 
ground state. Yet Eq. (C3l cannot be the exact ground 
state of the problem at this value of N since it is not 
analytic in 7. However, it can be the exact ground state 
in the limit iV — > 00. Thus we conjecture that Eq. (C3l 



is the ground state of the Hamiltonian in our problem 
in the limit N ^ 00. If so, it implies that our problem 
undergoes a quantum phase transition as a function of 7. 



Let us consider the time independent version of the 
model studied in this paper. It is given by the time in- 
dependent version of the Hamiltonian Eq. ^ 

H = -'^b^b+^S' + ^{b^S- + bS+). (CI) 



N 



In the sector where 



b^b + S' 



N 
2"' 



which is the one mostly studied in this paper, this Hamil- 
tonian can be thought of as simply aA^+1 byA^+1 
matrix. In the basis of states \n) where n denotes the 
boson occupation number, this matrix takes the form 



E^n'n 



9_ 
N 



9 



l\/N ~n' 5n' + l,n + 



nSn'-l, 



(C2) 



(compare with Eq. ([9|). 

It is easy to evaluate the eigenvalues of this matrix 
numerically for moderate N. 
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FIG. 16: The energy levels (eigennvalues of the matrix 



Eq. (C2|) for A'' = 10 are plotted as a function of 7. At the 
same time, the minimum of the classical Hamiltonian, given 
by Eq. ((C3|, is also plotted (red curve). We see that Eq. ([CS]) 



closely follows the ground state energy of Eq. ( C2 1 



When 7 depends on time, it follows that we drive our 
system across the quantum phase transition. The exis- 
tence of the critical trajectory and the singular behavior 
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of the adiabatic invariants discussed in section IVTl can be 
traced to the existence of this transition. The transition 
in the time independent version of the Dicke model is 
well known in the literature (see, for example, Ref. [46 
and references therein). 



At larger e, we can estimate it by 

Here, cq represents the moment in time when the sys- 
tem crosses the critical trajectory. At this point n did 



APPENDIX D: DERIVATION OF EQ. (62| 



Here, just as in the section VI we use the rescaled vari- 
ables n and H, according to Eqs. (39 1, (40 1. To obtain 
the value of e = e* at which w = tt, we need to inspect 
the increment of the angular variable w in the problem 
where 7 changes in time. It satisfies the equation [35] 



dtw ■ 



OA 
~dl 



where A was defined in Eq. ( 53 ) . At small A we neglect 
the second term to arrive at 



de' 

«o A log 



H 



(Dl) 



Here we used Eq. (58 1, the identity de/dt — 2A, and 
eo > denotes the moment of time when the system 
crosses the critical trajectory and switches to trajectories 
winding about the critical point Eq. ( 45 1 . This allows us 



to find w as a function of e, which in turn represents time. 
To compute this integral, we need to know H{e). This in 
turn can be estimated by noting that 

dtH = -2A71. 

Here, n itself is a function of time, which satisfies its 
equation of motion dtn 



which integrates to 



H 



-no 



, I 3A 

de e 



(D3) 



This needs to be substituted into Eq. (Dl I. We are now 



in a position to estimate the value e* corresponding to 
ui = TT. To compute the integrals in Eq. (Dll, we need 



to study the behavior of H. At small e — eg, H can be 
approximated as 



-no (e - eo) • 



not have a chance to change appreciably from its initial 
value, n{j — s- —00). Thus tiq ~ ^. This means that 
to logarithmic accuracy, we may use the approximation 
H{e) ^ —1/N in (Dl I. Other factors, such as 16e^ under 



the logarithm can also be neglected. This leads to 



de 



Alog[iV] 



2e*i 
3AlogiV 



Solution of the this equation for e* obtains the result 
( 62 1 . Our derivation has been self-consistent in the sense 



that even at the maximal value e* 



As 

;iiog^)- 



|iJ(e* 



where again the conditions Eq. ( |37| were used. 

The last thing which remains to be done is to check 
that Eq. (60) is consistent with the behavior of the sys- 
tem. The equation of motion Eq. (59 1, together with the 
condition that \H\ ^ e^, implies 

2Ade(/' = -e-h?!)^ 

assuming that (j)'^ < e. This is a Riccati equation, which 
can be brought to the form of the Schrodinger equation 
by the substitution 



2n(j). Using Eq. (|6()l we find 

(D2) to give 



R 



- f <pde 



d^R 

-2\-— + eR^Q. 
de^ 



This is the equation for the Airy function. It can be inves- 
tigated using the WKB approximation, which reproduces 
the ansatz (f) = —y/e. Close to e = 0, the WKB approxi- 
mation breaks down, so this ansatz is no longer correct. 
However, it is easy to modify relations such as Eq. (D2| 
and Eq. (D3l by substituting i?^, with R being the ap- 



propriate Airy function, in place of the exponentials in 
these relations. This does not affect the final answer for 
e*, and consequently for I final- 
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